Strongly correlated fermions after a quantum quench 
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Using the adaptive time-dependent density-matrix renormalization group method, we study the 
time evolution of strongly correlated spinless fermions on a one-dimensional lattice after a sudden 
change of the interaction strength. For certain parameter values, two different initial states (e.g., 
metallic and insulating), lead to observables which become indistinguishable after relaxation. We 
find that the resulting quasi-stationary state is non-thermal. This result holds for both integrable 
and non-integrable variants of the system. 
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Recent experiments on optical lattices have made it 
possible to investigate the behavior of strongly corre- 
lated quantum systems after they have been quenched. 
In these experiments, the system is prepared in an ini- 
tial state |"0o) an d then is pushed out of equilibrium by 
suddenly changing one of the parameters. Prominent ex- 
amples are the collapse and revival of a Bose-Einstein 
condensate (BEC) [1|, the realization of a quantum ver- 
sion of Newton's cradle I2j, and the quenching of a fer- 
romagnetic spinor BEC [3(. All of these systems can be 
considered to be closed, i.e., have no significant exchange 
of energy with a heat bath, so that energy is conserved 
to a very good approximation during the time evolution. 
Furthermore, since these systems are characterized by a 
large number of interacting degrees of freedom, appli- 
cation of the ergodic hypothesis leads to the expectation 
that the time average of observables should become equal 
to thermal averages after sufficiently long times. Various 
authors have recently given voice to such an expectation 
0, 0, 0]- However, the experiment on one dimensional 
(ID) interacting bosons shows no thermalization, a be- 
havior that was ascribed to integrability Q ■ Rigol et al. 
found that an integrable system of hard-core bosons re- 
laxes to a state well-described by a Gibbs ensemble that 
takes into account the full set of constants of motion Q; 
similar results were found for the integrable Luttinger 
model 0. 

For a closed system, the set of expectation values of 
all powers of the Hamiltonian H constitute an infinite 
number of constants of motion, irrespective of its inte- 
grability. Therefore, the question of the importance of 
integrability in a closed system that is quenched arises. 
We address this issue by investigating the full time evo- 
lution of a strongly correlated system whose integrability 
can be easily destroyed by turning on an additional in- 
teraction term. We show, using the recently developed 
adaptive time-dependent density matrix renormalization 
group method (t-DMRG) [|, DB EI El, that, in 
a certain parameter range, two different initial states 
with the same energy relax, to within numerical preci- 



sion, to states with indistinguishable momentum distri- 
bution functions. A comparison with quantum Monte 
Carlo (QMC) simulations shows, however, that they do 
not correspond to a thermal state. By using a generalized 
Gibbs ensemble @, @, HI, IH with the expectation value 
of the powers of the Hamiltonian (H n ) as constraints, we 
can improve the agreement with the time averages of the 
evolved system. This applies to both the integrable as 
well as the non-integrable case. 

In this Letter, we investigate the Hamiltonian 
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with nearest-neighbor hopping amplitude ih and nearest- 
neighbor interaction strength V at half-filling. The 
annihilate (create) fermions on lattice site i, n, = c\cj, 
and we take K = 1. We measure energies in units of 
ih, and, accordingly, time. The well-known ground- 
state phase diagram for the half-filled system consists 
of a Luttinger liquid (LL) for V < V c = 2, separated 
from a charge-density- wave (CDW) insulator (V > V c ) 
by a quantum critical point V c [161 ], This model is in- 
tegrable, with an exact solution via the Bethe ansatz 
[171 ]. We consider open chains of up to L = 100 sites 
pushed out of equilibrium by suddenly quenching the 
strength of V from an initial value V(t = 0) = Vq 
to a different value V(t > 0) = V. Furthermore, we 
study the effect of adding a next-nearest-neighbor repul- 
sion Vi ^2 ■ njTij + 2 to the model, which makes it non- 
integrable. We compute the time evolution using the 
Lanczos time-evolution method 13|, ID, lj|, 2(| and the 
adaptive t-DMRG. We study the momentum distribution 
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function (MDF) (n k ){t) = ± ££ r 

i.e., the Fourier transform of the one-particle density ma- 
). In the t-DMRG, we utilize the Trot- 



trix, pi r 
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tcr approach developed in Refs. [l(j as well as the 
Lanczos approach 12, E3l with additional intermediate 



time steps added within each time interval [21( . We hold 
the discarded weight fixed to e < 10~ 9 during the time 
evolution, but additionally restrict the number of states 
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kept to be in the range 100 < m < 1500. In all calcula- 
tions presented here, the maximum error in the energy, 
which is a constant of motion, is 1%, and, in most cases, 
less than 0.1%, at the largest times reached. 
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FIG. 1: (color online) Time evolution of an initial LL state 
(Vo = 0.5) in the strong coupling limit with V = 100 for 
L — 100 sites at the times indicated. Inset: spectral analysis 
of (n,r(t)}. 

When an initial LL state is quenched to the strong- 
coupling regime, V » t^, we find that (rifc(i)) (Fig. 
Q]), exhibits collapse and revival on short time scales, 
whereas the density-density correlation function remains 
essentially unchanged, i.e., retains the power-law de- 
cay of the LL. This can be understood by considering 
a quench to the atomic limit, th = 0. In this limit, 
all observables that commute with the density operator, 
including the density-density correlation function, arc 
time- independent. Furthermore, since the only remain- 
ing interaction is the nearest-neighbor density-density 
interaction, it can be shown analytically that the one- 
particle density matrix p m i (t) involves only two frequen- 
cies (loi = V and ui2 = 2V), resulting in a periodic oscil- 
lation with a revival time of T rcv ; va i = 2tt/V j22| . Thus, 
in analogy to the observed collapse and revival of a BEC 
in an optical lattice [l|, the single-particle properties of 
an initial LL state exhibit collapse and revival with this 
period. For the strong-coupling regime, the time evolu- 
tion retains the oscillatory behavior of the atomic limit; 
two frequencies u>\ and u>2 are indeed dominant in the 
spectrum, as can be seen in the inset of Fig. [1] However, 
the finite hopping amplitude leads to a dephasing of the 
oscillation on a time scale of idephase ~ 1/th- 

Afterwards, observables oscillate with a small ampli- 
tude around a fixed value, suggesting that the system 
reaches a quasi-stationary state. In order to further char- 
acterize such states, we study the evolution of the system 
for various values of V up to times one order of magni- 
tude larger than 1/th when applying the t-DMRG and up 
to two orders of magnitude larger when using full diago- 
nalization (FD). We find that the time averages for the 
longer times reachable by the FD agree with the time av- 
erages for the times reachable by the t-DMRG. Therefore, 
we conclude that the relevant time scale for the relaxation 
is indeed given by 1/th- In Fig. [21 the MDFs, obtained 
by performing an average in time from time t = 3 to 



t = 10 at the quantum critical point, V = V c , and at a 
point in the CDW region, V = 5, are shown. In order to 
investigate to what extent the (quasi-)stationary behav- 
ior is generic, we examine its dependence on the initial 
state. We do this by preparing two qualitatively different 
initial states with the same average energy (H) for each 
case: one a ground state in the LL regime and the other 
a ground state in the CDW regime. This is possible for 
a certain range of V in the intermediate coupling regime. 
In Fig. [2j results for two such initial states are compared 
with each other and with the MDF obtained for a sys- 
tem in thermal equilibrium and the same average energy, 
calculated using QMC simulations [23I ]. 

At the critical point, V = V c , Fig. G^a), the MDFs for 
the two initial states coincide with each other, to within 
the accuracy of the calculations (approximately the sym- 
bol size) or less. Therefore, information about the initial 
state is not preserved in this quantity, consistent with the 
expectation for an ergodic evolution. However, the dif- 
ference from the thermal distribution is significant; ther- 
malization is not attained. The left inset shows the time 
evolution of (rik}(t) for k = 7r for both initial conditions, 
demonstrating that the system reaches a quasi-stationary 
state. A small, but discernible shift from the thermal 
value (horizontal line) can also be seen for both initial 
conditions, even at k = n. In the right inset, the points 
with the largest differences just below the Fermi vector 
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FIG. 2: (color online) Time-averaged momentum distribu- 
tions when quenching (a) from Vo — 0.5 to V = 2 (quantum 
critical point) and (b) from Vo = 1-5 to V = 5 (insulator) for 
L = 50 sites. The time averages of two independent initial 
states with the same energy are compared to each other and 
to the thermal expectation value. In the right inset, results 
for L = 50 (*) are compared to L = 100 (v) for the regions 
with the largest differences. As a reference, finite T data for 
L = 50 (dotted line) and L = 100 (dashed line) are shown. 
Left insets: (n^) vs. time t; the horizontal line is the finite T 
value. 
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FIG. 3: (color online) Time- averaged momentum distribution 
(rife) for V=10. Inset: (n T ) vs. time t; the horizontal line is 
the finite-T value. 



kp = 7r/2 are plotted when the system size is doubled 
from L = 50 to L = 100. An analysis of the data for 
L = 100 does not alter our conclusion. 

In order to investigate the importance of quantum crit- 
icality, we have also examined the behavior for V = 1.5 
and V = 2.5 (not shown), i.e., below and above the tran- 
sition point. We find almost identical behavior, indicat- 
ing that the lack of thermalization is not associated with 
the quantum critical point. Note, however, that the LL 
regime (V < 2) is, in a sense, generically critical. 

For V = 5, Fig. [5] (b), all three curves show small, but 
significant differences. This means that the time evo- 
lution starting from different initial states can be dis- 
tinguished from each other as well as from the thermal 
state, i.e., neither relaxation to one distinguished quasi- 
stationary state nor thermalization occurs. For this case, 
| Vo — V | is larger than for V = 2, and, in addition, the 
values of Vq necessary to obtain the same energy in the 
initial state (Vo = 1.5 and Vo = 44.2165) differ strongly. 
This suggests that the initial states are far apart from 
each other in some sense, a notion that will be made 
more precise below. 

The differences with the thermal distribution increase 
for larger | Vo — V |. As can be seen in Fig. [3]for V = 10, 
the difference between the time average and the thermal 
distribution is significant, clearly confirming that ther- 
malization does not occur. The differences observed in- 
crease gradually as a function of | Vo — V |, ruling out a 
transition as suggested for the Bose- Hubbard model Q . 




FIG. 4: (color online) Comparison of QMC and time-averaged 
values of the momentum distribution (rik) for the non- 
integrable case with V — 2 and Vi = 0.4. Inset: the same 
as in Fig. [3] 

In order to investigate the impact of the lack of in- 



tcgrability on thermalization, we now extend our model 
(fTj) by turning on a next-ncarcst-ncighbor interaction. In 
Fig.gl we display results with Vo = 0.5 and V = 2.46689 
(zero n.n.n. interaction), and the quenched evolution at 
V = 2, Vz = 0.4. As in the integrable case, both ini- 
tial states lead to indistinguishable time-averaged MDFs, 
but ones that are significantly different from the thermal 
one, showing differences very similar to those in Fig.[2ja). 
When V = 0.5 and V = 10, V 2 = 1 (not shown) the dif- 
ference from the thermal state is comparable to the case 
shown in Fig. [3] Therefore, the non-thermal nature of 
the emerging steady state is clearly not related to the 
integrability of the system. 

In order to shed light on the numerical results pre- 
sented above and to characterize the quasi-stationary 
state, wc consider a generalized ensemble in which the 
expectation values of higher moments of H, which are 
constants of the motion, are taken as constraints. Note 
that the usual thermal density matrix gp is uniquely fixed 
by the single constraint (H)p = (H). Rigol et al. 0], 
find a generalized Gibbs ensemble, in which the density 
matrix is determined by maximizing the entropy taking 
into account the constraints, to be an appropriate choice 
0,[l4, 15 1 . The general form of the statistical operator 

is then g = exp — J2 n ^nO n , where the operators O n 
form a set of observables whose expectation values re- 
main constant in time. The values of the X n are fixed by 

the condition that Tr [gO r ^j = (O n ), with Oq = 1 to en- 
force normalization. In some special cases like hard-core 
bosons in one dimension 0, HH or the Luttinger model 
[H , constants of motion can be found in terms of opera- 
tors in second quantization. However, this is not possible 
for Bethe-ansatz-integrable systems. For any closed sys- 
tem, however, the quantities O n = H n can be used. Tak- 
ing all powers as constraints would unambiguously fix all 
correlation functions to all lengths. For a finite system, it 
can be shown that g is fully determined by dim(if) pow- 
ers of H |22j, for H with a bounded spectrum. The sta- 
tistical expectation value of any observable is then given 

by Tr [gOj = X) l( i/ l^'o)| 2 (HO|i') (for a non-degenerate 

spectrum), where |"0o) is the initial state and \v) are the 
eigenstates of H. It can be easily seen, that the r.h.s. of 
this expression equals the time average of (O) (t) . 

We now investigate the extent to which the statistical 
expectation value within the generalized Gibbs ensem- 
ble approaches the time average of the evolution after a 
quench by studying the energy distribution for a given 
state defined as P 4 ,{E) = S (E - t v ) \ (v\ip) | 2 , 
which is normalized if (V'lV') = 1- The energy distribu- 
tion in the generalized Gibbs ensemble can analogously 

be defined as Pq (E) = Tr 5 (^E — H^j g, with g as defined 

previously. In Fig. El we show P^(E) calculated using 
full diagonalization on L — 16 sites for an initial state 
Vo = 0.5 evolved at the quantum critical point, V = 2 
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V„ = 0.5.V = 5 




lost in the observables after relaxation. 'Closeness' can 
be quantified using a measure Ao which is based on the 
energy distributions defined for the initial state or for 
a given density matrix. Increasing the number of con- 
straints (moments of H) in a generalized Gibbs ensemble 
leads to convergence to the energy distribution defined 
by the initial state. 



FIG. 5: (color online) Influence of the constraints (H n ) on 
the energy distribution function P(E). 



compared with the distribution in the Gibbs ensemble 
Pg(E) as the number of constraints is increased from 1 
to 3. It is evident that increasing the number of con- 
straints systematically improves the agreement and that 
only a small number of moments are necessary to obtain 
very good agreement. 

The distance between two distributions can be esti- 
mated using the moments of the absolute differences, 
A„ = / dEE n | P(E)-P'{E) | . Taking A /W, with W 
the bandwidth of H, as an estimate of | P(E) — P'{E) |, 
we see that the difference between moments of H for two 
different energy distributions P and P' can be estimated 
as 



(H n ) P - (H n ) P , < A n 
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Therefore, if the distance between the distributions 
Ao C 1, then the relative difference of the moments 
({H n ) P - (H n ) P ^ /W n < A /(n + 1) will also remain 

small, and observables will converge to values close to 
each other after a quench. For the cases of evolution with 
metallic and insulating initial states discussed above, we 
obtain A = 0.12439 for V = 0.5 and V = 3.57463 (V = 
2), and A = 0.41521 for V = 1.5 and V = 44.2165 
(V = 5). On the other hand, comparison of P^ with the 
thermal distribution Pp yields A = 0.68581 (V = 0.5, 
V = 2), and A = 1.24616 (V = 1.5, V = 5), respec- 
tively. Thus, the distance between the thermal distri- 
bution and the one defined by the initial states is always 
larger than those defined by the pair of initial states with 
the same energy, supporting our observation that ther- 
malization does not occur. 

In summary, our adaptive time-dependent density- 
matrix renormalization group simulations of the time 
evolution of a system of correlated spinless fermions after 
a quantum quench have exhibited the following generic 
behavior: Independently of its integrability or critical- 
ity, the system relaxes to a non-thermal quasi-stationary 
state. Observables relax to the same value when differ- 
ent initial states have the same energy and are sufficiently 
close to each other, i.e., the memory of the initial state is 
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